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Abstract 

We explore the praticability of optimal shape design for flows modeled by the 
Euler equations. We define a functional whose minimum represents the optimality 
condition. The gradient of the functional with respect to the geometry is calculated 
with the Lagrange multipliers, which are determined by solving a costate equation. 
The optimization problem is then examined by comparing the performance of several 
gradient-based optimization algorithms. In this formulation, the flow field can be 
computed to an arbitrary order of accuracy. Finally, some results for internal flows 
with embedded shocks are presented, including a case for which the solution to the 
inverse problem does not belong to the design space. 
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1 Introduction 


A classical problem in engineering is to define the shape of a vehicle to achieve a required 
performance. In fluid dynamics, techniques have been developed to solve the following 
inverse problem: given a pressure or a velocity distribution over an aerodynamic body, 
determine the corresponding geometry. See, for example, reference [8]. A broader category 
of problems can be solved by means of optimization, provided that one is ready to accept 
the necessity of computing the flow field hundreds of times. In using models of increased 
complexity to describe the flow field, the development of new algorithms is necessary in 
many cases to reduce the computational load. In this paper, we investigate one method for 
achieving this reduction. The variational technique that is applied in this work has been 
used since before complex flows could be integrated numerically. See, for example, refer- 
ence [1]. Jameson [4] was the first to apply this technique to computational fluid dynamics. 
With this approach, a functional or cost function is determined such that its minimum 
represents an optimal solution. By introducing a set of Lagrange multipliers, the gradient 
of the functional can be calculated with respect to the geometry by computing the flow 
field only once for each gradient evaluation. For incompressible irrotational steady flows, a 
further reduction in the computational effort is possible. (See [10].) The formulation devel- 
oped in this work is more suited to a correct analysis. In fact, this formulation eliminates 
the need to consider the flow field variables dependent on the geometry. In the present 
work, we extend the work presented in reference [5]. In reference [5], an exact gradient 
of the functional with respect to the design variables was obtained on the discrete level, 
which can be a limitation for compressible flows because in presence of shocks the discrete 
functional can present discontinuities. In considering compressible flows with embedded 
shocks, we derive the gradient on a continuous level. Furthermore, we provide a method 
for calculating the conditions that the Lagrange multipliers must satisfy at the boundaries 
and at the shock. Finally, we point out that our formulation can be used with complex 
flow solvers because the differentiability of the solver is not requested. 


2 Problem Statement 

The Euler equations are given by 

U t + F x + G y = 0 (1) 

where 
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u = x component of velocity vector 
v — y component of velocity vector 
e = specific total energy 
pressure 
speed of sound 
specific heats ratio 
7 ~ 1 
2 

and p = np( 2e — u 2 — v 2 ). Furthermore, 

F=|£u = A(U)U (2) 


P = 
a = 

7 = 

K = 


and 


/9C 1 

G = -U = B (U) U 


(3) 


where A and B are given in app. I. 

We assume that these equations are defined on a physical space $. In this space is included 
a subdomain f! whose boundary is denoted by T. On the boundary, we define a curvilinear 
coordinate s and a normal n = (n x ,n y ) that points outward. (See fig. 1.) 

The optimization problem studied here is defined as the minimization of the functional 
£ = f r (f>(p, p,u,v,r)ds over all admissible shapes of the subdomain fl, subject to the 
steady-state Euler equations with proper boundary conditions on T. 

Although the method we present is general, we focus on the following model problem. The 
subdomain 0 is represented by a nozzle. (See fig. 2). At the inlet, total pressure, total 
temperature, and the ratio a = vju are fixed. At the outlet, if the flow is subsonic, the 
static pressure is fixed, and at the solid walls the impermeability condition un x + vn y = 
0 is enforced. The upper wall is kept fixed. The lower wall 0 is represented by the 
parameterization 

J(6) = l>/i(*) (4) 


where the functions fi{x) are some shape functions and a = (aj, ..., a,-, ...) is the corre- 
sponding set of shape coefficients. Given a desirable lower wall pressure distribution p*(x) 
and the actual pressure distribution on the lower wall p u '(x), the optimization problem 
consists in finding a set of shape coefficients a,- such that the functional 


£ = 



(5) 


is minimized. 


2 



3 Lagrange Multipliers and Optimality 

The problem of achieving the minimum is addressed by introducing a set of Lagrange 
multipliers. Consider the augmented functional 

£(U, a, A, //) = £ + f t A( AU X + BU y )dfi + f ppV • n ds (6) 

JQ J 0 

where V = ( u , u). The vector A(x, y ) = 4 (Aj, A 2 , A 3 , A 4 ) and the scalar p(s) are the contin- 
uous equivalents of the Lagrange multipliers. 

We calculate the variation of the functional C with respect to the variation of the functions 
U, A, and p and the parameters cq , respectively. When U(x,y) is increased by a function 
eU (x,y), the functional C increases by an amount eSCjj. In the same way, A(x,y) is 
increased by eA (x,y); p(s), by ep(s)-, and each a,, by ea t . 

If we follow the derivation in app. II and take 


SC — SCu d* SC\ T SC.fi 4- SC a 


then we obtain 


SCu = I 6 ~ ( p w - p')Vdz + [ *A(A n x + Bn y )Uds + 

Ja U U JT 


- f ( l A x A + 'A y B)Udfl + [ p. n ^ U ds 
Jn J@ a U 
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Furthermore, 

SC K = / ‘A(AU* + BU y )<m 

Jq 

SCu = / ppV ■ nds 

J 0 

SC a = y\Ft- W U -P*)fidx+ / ( A(AU X + BU y ) /, cos 9 ds 

i |/a dy q JQ 


+ / P 


f p - -nfids- f ppV ■ t 
Je oy Je 


© dy 


dx 


cos 2 6 ds 


+ / /ipV ■ n — - sm9 cos 9 ds d; 

Je dx 

where 9 is the angle between the normal n and the y-axis and t = (— n y ,n x ). 


(7) 


( 8 ) 

(9) 


( 10 ) 
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At the minimum of the functional, we have SC = 0 for all possible choices of the functions 
U, A, and p and of the parameters a. This condition is reached when 

SCu = SC\ — SC M = 6C a = 0 (11) 

Note that because of the necessary conditions (eqs. (11)) the unconstrained minimum of the 
functional £(U,a, A,/r) corresponds to the constrained minimum of the functional £{a). 
In fact, we have 

SC\ = 0 => AU X T BU y = 0 and SC ^ = 0 => p\ ■ n = 0 

which means that U must satisfy the Euler equations with boundary conditions. In addi- 
tion, for the minimum of C we have SCu — 0, which leads to 

*AA X + ‘BA,, = 0 on ft (12) 

= 0 on 0. (13) 

At the inlet, outlet, and upper wall, 

‘A(A n x + Bn y )U = 0 (14) 

Given U and the set of costate eqs. (12) and (14), we can determine uniquely A in and 
p on 0. (See app. III). 

Finally, given a and given U and A from the above equations, we can calculate from eq. 
(10) 

(»* -»•)/( ir + f 'A(AU. + BU,)/, cos (Ids 

+ f (U— ^ — - • n fids — f ppY ■ t -p- cos 2 0 ds 
Je ay Je dx 

+ f pp\ ■ n sin 0 cos 0 ds (15) 

Je dx 

In cases for which shock occurs in the flow field, we split the domain of integration by 
means of a curve T that coincides with the shock where it exists. Then, we follow the same 
derivation so far on each of the two subdomains, with T as a boundary. (See app. III.) 
The strategy that we use to achieve the minimum of C is as follows: 

1. Start with a set a of shape coefficients. 

2. Enforce SC\ = 0 and <5£ M = 0 by finding U such that it satisfies the steady-state 
Euler equations and boundary conditions. 

3. Enforce SCu — 0 by finding A such that it satisfies the costate equations and 
boundary conditions. 


dC _ r b dp 
da i J a dy 0 


d£_ /_« 

T \P 


p*) cos 6 + ( A( A n x + Bn y ) + /m 
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4. Calculate V Q £. If V a C = 0 then we have determined the minimum; otherwise 
continue to steps 5 and 6. 

5. Update a with criteria based on V Q £. 

6. Restart from step 2. 

4 Discrete Problem 

We introduce a discrete grid that is defined as ( x\,y m ) = ( x 0 + /Ax,y(0) + mAy), where 
Ax is constant and Ay is a constant fraction of the local height of the nozzle. (See fig. 3.) 
The steady solution of the Euler equations is obtained with a time-dependent technique, 
in the frame of an explicit finite- volume code. The conservative variables U are computed 
at the cell centers, and the fluxes F and G are evaluated at the cell interfaces with the 
approximate Riemann solver in reference [7]. Second-order accuracy is achieved by using 
an essentially nonoscillatory scheme [2]. The flow-field values at the cell interfaces, used 
as initial conditions for the Riemann problem, are reconstructed by means of a linear 
interpolation and a minmod limiter. The amplitude of the integration step is chosen in 
accordance with the Courant-Friedrichs-Lewy (CFL) condition. 

The costate equations are discretized on the same grid presented above. Because they have 
no conservative form, the numerical solution is obtained with a finite-difference scheme. 
We introduce a set of curvilinear coordinates tp(x,y) and rfr(x,y). The costate equations 
are then written 

t AA v + t BA 4> = 0 (16) 

where A = A <p x + Bc^ and B — A xj> x + B^v The transformations <p and V 7 are defined as 
(x h y m ) / and {x h y m ) ^ m, respectively. 

We find the solution of eq.(16) as the asymptotic limit of a time- dependent technique. 
Consider eq. (16) embedded in time as 

± A t + * AA V + t BK 4 , = 0 (17) 

We must select the proper sign for the time derivative. The inlet and outlet boundary 
conditions for the costate equations are complementary to those of the flow field equations, 
in the sense that if the number of boundary conditions for the flow field is c, then the number 
of boundary conditions for the costate equations is 4 — c. Therefore, the above equations and 
boundary conditions are well posed if we select the negative sign for the time derivative. 
In fact, the resulting characteristic pattern is mirror symmetric with respect to that of the 
flow- field equations. 

In the presence of a shock in the flow field, the matrices l A and *B are discontinuous. In 
particular, the characteristic pattern at the shock indicates the necessity of a boundary 
condition for the costate equations. For further discussion, see reference [3]. 

The costate equations are linear and as such are the boundary conditions. We exploit 
this property to numerically solve these equations. Suppose that locally we separate the 
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( 18 ) 


variables through the following approximation: 

A (v?,tM) = A'(vM) + A "(tM)- 


This separation of variables means that, for example, in a Taylor expansion about the point 
(<p,rp) we disregard all terms that involve a cross product iptp. This approximation is at 
least first-order accurate. We substitute eq. (18) into eq. (17) to obtain 

- a' ( -a" ( +x + ( m; = o 

and we are left with the following subproblems in one dimension: 

-A\+Ua; = 0 (19) 


- A" t + l BAl = 0 (20) 

Let us define n v = (Vx/fal + vliVy/ \Jvl + ¥>*) and = (tp x / yty* + ip y / f + ^). 

The left and right eigenvector matrices of A and B are calculated by using the formulas in 
App. I with n = and n = n^, respectively. After eqs. (19) and (20) are diagonalized, 
we upwind the derivatives of the characteristic variables according to the signs of the 
corresponding eigenvalues. The time step At is chosen according to the CFL condition. 
This method can be regarded as a two-dimensional interpretation of the method presented 
in reference [6]. 

The boundary conditions can be split in a similar way. Consider, for example, the boundary 
condition at the solid wall. Because eq. (19) is defined along the wall, the characteristic 
variables can be upwinded according to the corresponding eigenvalues. On the contrary, 
the third row of eq. (20), which corresponds to the characteristic with a speed of +a, is 
replaced by the boundary condition in eq. (AIII.5). Note that the contravariant component 
of the speed in the direction xp is 0; therefore, the resulting system can be written as 


att; 

A f W'' 

i n x A^Aj + n y A 4 A 3 

att; 


0 

0 

— (p w — p*) cos 0 — n x A^j — n y A t \' 3 

aA((v7|T^)A»< ' 


( 21 ) 


where A(-) is the forward finite increment of the function (-) with 
scripted variable and 


AW = 


AWi \ 

AW 2 

aw 3 

A W 4 / 


‘L^AA 


respect to the super- 


In the third row of eq. (21), we have the functions of A^, which are computed separately 
as mentioned. The other boundary conditions are enforced in the same pattern that is 
presented above. 
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5 Optimization Experiments 

The optimization problem is addressed with four different gradient- based criteria. 

1. Steepest descent (SDl). The shape coefficients are updated as follows: 
cti <— a, — v{dC/da t ). where v is a given parameter. 

2. Steepest descent with v selected as shown (SD2). Because we know the gradient 
V Q £ at the present iteration, we can use a tentative v and can compute the 
gradient V Q C ' . By calculating V a C ■ V a £', we linearly estimate v such that 
eventually V a £ • V 0 £” = 0. Each step of the optimization requires solution of 
the flow-field and costate equations twice. 

3. The BFGS algorithm (presented in reference [9]). This algorithm (BFGS1) ac- 
counts for the curvature of the hypersurface £. The shape coefficients are up- 
dated according to the formula a, *— c*i — where d = (...,d,,...) is the 
descent direction determined by d = H V a £ and H is an estimate of the inverse 
of the Hessian of the functional. 

4. The BFGS algorithm (as above) with a linear estimate of v such that d • V Q £" = 
0. Each optimization step requires solution of the flow-field and costate equations 
twice (BFGS2). 


The computations are performed on a 40 x 20 grid unless otherwise specified. Total pressure 
and total temperature at the inlet are taken unitary and cr(0, y) = 0. At the outlet, the 
static pressure depends on the test case considered. For the lower wall ordinate ?/(©), we 


have 


0 

y(@)={ £i=i c*iX ,+1 (x - l ) 2 
0 


if —0.5 < x < 0 
if 0 < x < 1 
if 1 < x < 1.5 


The optimization consists in finding the four shape coefficients a, such that the modulus of 
the gradient V a £ is 0. The flow-field and costate equations are iterated until the residuals 
are less than 10 -5 . 

Consider a test case in which a pressure distribution in a subsonic nozzle is found for which 
the outlet pressure is 0.9 of the inlet total pressure. We take a = (2, 0,0,0) and define the 
corresponding configuration as the optimal configuration. Then, we compute the flow field 
and determine the pressure distribution on the lower wall. This pressure distribution p* is 
the one we want to recover with the optimization algorithm. In fig. 6, the target flow field 
and the starting configuration, obtained with a = (—2,0, 0,0), are shown along with the 
convergence histories for DS1 and BFGS2. 

For the supersonic case, we take a constant section channel as a starting configuration and 
a — (2,0, 0,0) as the target. In fig. 7, we present the results obtained when the outlet 
pressure is lowered to 0.5 of the inlet total pressure. A relevant shock is generated in the 
flow field. In the first optimization iterations, we updated the shape coefficients, as was 
done for DSl. This step is necessary because this far from the minimum the functional 


7 



£ might be not convex; therefore, the estimate of v that was used would not be correct. 
Figure 8 shows the sequence of lower wall configurations obtained with BFGS2. We also 
tested the effect of grid coarseness by reducing the number of grid points to 20 x 10. (See 
fig- 9.) 

The second test case is designed to check the capability of the algorithm in detecting 
minima in cases for which the desired pressure distribution is out of the design space (i.e., 
the functional is not 0 at the minimum). The pressure distribution p* is obtained with an 
outlet boundary condition that differs from the distribution that is actually used in the 
optimization routine. The results are given in fig. 10. 

The DS1 updating strategy had the least attractive rate of reduction of the functional. Our 
experience showed, nevertheless, that it was the most reliable in cases of complicated surface 
topologies that can occur in flow fields with embedded shocks. The BFGS2 becomes the 
most efficient of the algorithms tested when it is coupled with DSl. With this algorithm, 
the functional was reduced by orders of magnitude. 

Each optimization procedure with a 40 x 20 grid required 6 h of central processing unit 
(cpu) time on a DEC 3000/500. The 20 x 10 case required 1 h of cpu time. 


6 Concluding Remarks 

We have derived an expression of the gradient of the cost function with respect to the shape 
coefficients. The boundary conditions for the costate equations have been presented; we 
have shown their relevance to the well posedness of the problem. In the case of shocks, 
we provided the proper conditions for the costate equations at the discontinuity. On the 
discrete level, we proposed a method of integrating the costate equations in accordance with 
a revisited scheme. Additional work is needed to test the algorithm with more realistic test 
cases and to apply the One-Shot method (reference [10]) to hyperbolic problems. 
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Appendix I 

The Jacobian matrices for the Euler equation in conservative variables are 



0 

1 

0 

0 

kV 2 - u 2 

(3 - 7 )u 

—2kv 

2k 

—uv 

V 

u 

0 

(—'ye + 2 kV 2 )u 

7 e — kV 2 — 2 ku 2 —2 kuv 

7U 

0 

0 

1 

0 

—uv 

V 

u 

0 

kV 2 - V 2 

—2ku 

(3 - 7)v 

2k 

(— 'ye + 2 kV 2 )v 

—2kuv 'ye — kV 2 — 2 kv 2 

■yv 


(AI.1) 

(AI.2) 


The Jacobian in the direction n is C = An^+Briy. The left (L n ) and right (I^ 1 ) eigenvector 
matrices of C are 


1 - kV 2 /a 2 2 ku/a 2 2 kv/a 2 -2k /a 2 ' 

Vt/P n y/p ~n x /p 0 / AT ,x 

( — V n + kV 2 /a)/p (n T — 2ku/a)/p ( n y — 2kvja)jp 2kjpa 
(V n + kV 2 /a)/p -(n x -\-2kuja)jp — (n y + 2kv/a)/p 2kjpa 

1 0 pj2a 

u pn y p(u + an x )/2a p( 

v pn T p(v + an y )/2a p( 

V 2 j 2 -pV t p(V 2 + a 2 /k + 2aV n )/4a p(V 2 + 

where V n = V • n and V t = V • t. The diagonal matrix D„ = 

V n 0 0 0 

ok o o 

0 0 V n + a 0 

0 0 0 V n — a 

Appendix II 

To calculate 6Ca , consider the increment U <— U + eU =$• 

F <— F + eAU + h.o.t. and G <— G + eBU + h.o.t. 

We obtain 

SC V = ^ (p - p*)Udx + jf ‘A [(AU) e + BU)„] dSl 

+ / pn Uds + h.o.t. ( A 1 1 . 1 ) 
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If we apply Gauss’ theorem to the second integral of the above equation then we find eq. 
(7). Equations (8) and (9) are easily obtained. 

To calculate 6£ a , we first consider the variation of the functions defined on 0: 

f O \ T 

ai <— a t + £&i =» p w *- p w + e-j- fi &i and pV <— pV + e fi cq. 

dy dy 

Then, we consider the variation of the geometry; other higher order effects are disregarded. 
When the geometry is perturbed, the domain of integration ft, the normal n and the 
element of integration ds are perturbed. The domain ft is increased (fig. 4) by a quan- 
tity edif l cos 8 ds. The normal n is perturbed by a quantity —eaidfi/dx cos 2 0t; ds, by 
eaidfi/dx cos 8 sin# ds. (See fig. 5.) 

Appendix III 

Consider eq. (14). This equation defines the boundary conditions for A after we impose 
the proper constraints on U. At the inlet, only one component of the variation of the flux 
in the direction normal to the boundary F„ = (An* + Bn y )U is independent of the others 
because total pressure, total temperature, and cr are fixed. If we express all components of 
F n in terms of pu, we obtain 

* 

[i + w/( 1 - M 2 )] u 

- vlh - m 2 )]u 

Hi 

where £ = n x + an y , p = (n y — an r ), H is the specific total enthalpy, and M is the local 
Mach number. For an arbitrary choice of pu, from eq. (14) we have 

Ai£ + A 2 u [£ + 77<r/(l - M 2 )] + A 3 u [<r£ - ??/( 1 - M 2 )] + A 4 /f£ = 0 (AIII.l) 

At the outlet, if the flow is subsonic, then only the static pressure is fixed; therefore, three 
components of the vector U are arbitrary. If we take pe as the dependent variable, we have 

pun x + pvn y 

pu(V n + un x ) + uriypv — puV n 
pv{V n + vriy ) + vn x pu - pvV n 
p{~;'p/2kp + u 2 + v 2 ) + pu(V n u -f Hn x ) + pv(V n v + Hn y ) 

For an arbitrary choice of p, pu, and pv, from eq. (14) we have 

A in* + \ 2 {V n + un x ) + A 3 vn x + A 4 (V n u + Hn x ) = 0 

Ai n y •+- \ 2 un y + As( V(j + vn y ) + \ 4 (V n v + Hn y ) = 0 (AIII.2) 

A2nV n + A3UI4 + A 4 V n ( r Yp/2kp + u 2 + v 2 ) = 0 
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For a supersonic outlet, no conditions exist on U; therefore, we have 


A = 0. (AIII.3) 

If a shock is embedded in the flow field, then the shock is considered as a boundary for 
the costate equations. The consequent boundary conditions are applied on each side of eq. 
(AIII.3) before the shock and eq. on each side of (AIII.l) after the shock. 

For the upper wall, we have 

0 ‘ 

P 

a y 

0 . 

such that for an arbitrary choice of p eq. (14) is satisfied if 



A 2 fi x + A 3 n y — 0 (AIII.4) 

At the lower wall, eq. (13) applies. Because the rank of A n x + Bn y at the wall is 2, the 
system has only two linearly independent rows. We obtain 


A 2 n x + A 3 n y + ( p w — p*) cosd = 0 
which is the boundary condition for A, and 


f i = — 


Ai + 11X2 + uA 3 + ( 7 e — fcV 2 )A 4 j 


which is the relation between p and A on the boundary. 


(AIII.5) 


(AIII.6) 
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Figure 1: Physical space. 



Figure 2: Model problem. 
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Figure 3: Discrete grid. 



Figure 4: Variation of domain of integration D with AH = dx , AD = ey and BC = 
e(y + dy/dx dx). 
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(a) Target Mach-number field. 


(b) Starting configuration. 



Functional 

igi 



Functional 

igi 


(c) Functional and modulus of gradient 
versus number of iterations for SDl. 


(d) Functional and modulus of gradient 
versus number of iterations for BFGS2. 


Figure 6 
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Figure 8: Wall shapes sequence with BFGS2. Starting configuration: constant section. 
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